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ABSTRACT 


The horizontal separations between convective thermal plumes, and 
features within a thermal plume, were determined through the use of an acoustic 
sounder, an anemometer and extensive data analysis. The mean, standard 
deviation, median and mode were calculated for the computed correlation 
lengths of the acoustic sounder data sampled in time intervals of 2, 5 and 10 


minutes. 


The data sampled at 2 and 5-minute intervals emphasized features within 
an individual thermal plume. The mean correlation distances found for 2 and 5- 
minute intervals was 81 meters + 70 meters and 89 meters + 72 meters, 
respectively. Their medians were 61 meters and 69 meters; and their modes 
were 41 meters and 50 meters, respectively. 


The 10-minute time interval statistics used a low pass filter to emphasize 
larger scale features. The mean correlation length was 494 meters + 373 
meters, the median was 391 meters and the mode was 316 meters. These 
distances represent the distance between the center of a plume and the center of 
a quiet region adjacent to that plume. 
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I. INTRODUCTION 


A. OBJECTIVE 

Light propagation through the atmosphere is affected by atmospheric 
absorption, scattering and turbulence. Velocity fluctuations in the atmosphere 
mix regions of the air with different temperatures producing local changes in the 
refractive index. This optical turbulence induces phase and amplitude 
fluctuations in electromagnetic waves that propagate through the air. These 
fluctuations smear the ability to see detail in imaging systems or spread and 
induce scintillation in laser beams. The atmospheric optical fluctuations are not 
distributed uniformly in the atmosphere, but are concentrated in local convective 
regions during the day. Understanding the spatial scale of these regions is part 


of the process of mitigating their effects. 


This thesis focused on measuring an atmospheric outer scale length 
known as the horizontal convective scale by computing the autocorrelation length 
of the acoustic sounder and anemometer data in order to estimate the lateral 
dimensions of the turbulent, detrimental events. The data analyzed in this 
thesis were generated by an acoustic sounder that was directed vertically into the 
atmosphere to collect turoulent atmospheric data and by an anemometer used to 
measure wind speed. These instruments were located at the Starfire Optical 


Range, New Mexico. 


B. BACKGROUND 

1. The Acoustic Sounder 

An acoustic sounder is a device designed to emit acoustic radiation and 
detect returns of that energy. This device has also been called a sodar (sonic 
detection and ranging). It is much like sonar, which detects objects underwater, 
except that sodar detects and measures properties of the air. One of the primary 
things that an acoustic sounder is good at is detecting and characterizing 


turbulence in the air since the signal that the sounder emits is scattered back by 


turbulence. 


Acoustic sounders are available commercially. An example of a 
commercial acoustic sounder is the Model VT-1 sodar developed by Atmospheric 
Research & Technology (ART) (Neal 2004). This unit is a self-contained, 
portable system and includes a phased-array acoustic transmitter and receiver 
with supporting electronics, a notebook computer and software for system 


configuration, operation and data storage. Figure 1 shows this sodar. 





Figure 1. Model VT-1 phased-array doppler sodar system. 


The acoustic sounder used to gather the data presented in this thesis was 
designed and built by Professor Donald Walters, Ph.D., professor of physics at 
the Naval Postgraduate School, Monterey, CA. Figure 2 shows a picture of this 
acoustic sounder. Acoustic sounders are not good at retrieving reliable data 
during rain and their internal components can be damaged by inclement weather. 


Notice the lid on the acoustic sounder which can be lowered during rain. 





Figure 2. Acoustic sounder designed and built by Professor Donald Walters 
and Jay Adeff at the Naval Postgraduate School. 


2. Convective Turbulence 
Thermal convective cells, or plumes, are formed in the atmospheric 
boundary layer (the layer of air bounded by the surface of the earth and about 1- 


2 kilometers above the surface) during the day as the sun warms the earth. 
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These convective cells consist of heated, turbulent air and they concentrate the 
optical path degradation into near vertical columns. These plumes are, however, 
separated by regions of relatively low turbulence where the optical path is 
relatively quiescent. Measuring the size of these plumes of turbulent air, and the 
quiet regions between them, is a component of an atmospheric compensation 
design. In order to understand and use an acoustic sounder for turbulence 
measurements, the structure of turbulence and the interaction with optical and 


acoustic radiation is essential. 


The following section, from equations 1 to 10, follows that done by 
Lim (2003). 


To deal with the randomness of atmospheric turbulence, 
Kolmogorov used a statistical approach that relies on dimensional 
analysis to handle the spatial and temporal fluctuations (Max 2003). 
By assuming homogeneity and isotropy at least in a local volume, 
and if the random processes have slowly varying means, structure 
functions represent the intensity of the fluctuations of f(r1,r2) over a 
distance between r; and re. Using the mean square differences, the 
structure function of f(r4,r2) is: 


D (rate) =([ 107) 102) ). (1) 


According to Kolmogorov’s turbulence theory, turbulent eddies 
range in size from macroscale to microscale, forming a continuum 
of decreasing eddy sizes. Energy from convection and wind shear 
is first added to the system at the outer scale L, (10’s - 100’s of 


meters) before it cascades to a smaller scale /, (~ 1cm) where 


viscosity converts the energy to heat (Andrews 2001). By 
dimensional arguments and assuming an incompressible isotropic, 
homogeneous medium, Kolmogorov showed that the longitudinal 
structure function of the velocity is: 


Disc, ao @ oe (2) 


The r’’® proportionality of the structure function in the inertial range 
(|, <r <L,) applies to other structure functions such as temperature 
and refractive index. The refractive index structure function is: 


D,(r) = C2?" (3) 


where C* is the refractive index structure parameter, or the optical 


turbulence parameter (m**%). Comparing equations (1) and (3), the 
refractive index structure parameter has the form: 


((a-ne)") 

5 oe 

While C? is the critical parameter that describes optical turbulence, 
it is extremely difficult to measure C* directly using standard 
techniques since the index of refraction of the atmosphere is 
influenced by the atmosphere’s temperature, pressure, moisture 
and the wavelength of the electromagnetic wave. However, C? 
depends on the temperature structure parameter, C* which can be 


measured directly. C? has a mathematical form similar to C*: 


ee, Cat) é 


2 
2/3 


Though the refractive index depends on the dry-air wavelength, it is 
typical to ignore the wavelength dependence and assume a 
wavelength of 0.5um (Beland 1996). The index of refraction is 


n=1+79x10°P/T. (6) 


Taking the partial derivative of the air density with respect to the 
temperature and assuming isobaric density fluctuations, the optical 
turbulence parameter C* relates to the temperature structure 


parameter C*> by 
on) PY 
o: - (| of -(7ox10* Z, | ox (7) 


where C? is the optical turbulence parameter in m®*, C? is the 


temperature structure parameter in K*m®%, P is the air pressure in 
mbar and T is the air temperature in Kelvin. 


3. Acoustic Sounder Measurement 


The acoustic sounder uses acoustic waves scattered by 
temperature and velocity fluctuations to measure changes in the 
refractive index of the atmosphere. It transmits an acoustic signal 
into the atmosphere and, for energy backscattered at 180 degrees, 


detects variations in the thermal structure parameter C? (Tatarski 
1971). Once C? is found, equation (7) provides the optical 
turbulence parameter C*. 


The power returned from the atmosphere was summarized by Neff 


(1975) 
ere} Ze)) ‘ 


where = is the received power (P. is the measured electrical 


r 


power and E, is the efficiency of conversion from received acoustic 
power). PE, is the transmitted power (P, is the electrical power 
applied to the transducer, E, is the efficiency of conversion to 


radiated acoustic power). eis the round trip loss of power 


resulting from attenuation by air where @ is the average extinction 
coefficient (m') to the scattering volume at range R(m). o, is the 


backscatter cross section per unit volume and S is the maximum 

effective scattering volume thickness where c is the local speed of 

sound (ms!) and z is the acoustic pulse length (s). eae is the 
R 


solid angle subtended by the antenna aperture A (m°*) at range R 
(m) from the scattering volume, modified by the factor G that 
accounts for the non-uniform antenna illumination. 


Tatarskii (1971) expressed the backscatter cross section o, at 180° 
as 





(9) 


where k=2z/A is the incident wavenumber, 7, is the mean 
temperature, and ®,(2k) is the 3-D spectrum of turbulence. The 
cross section o, represents the in-phase addition of backscattered 
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waves from temperature inhomogeneities spaced 1/2 apart along 
the radial propagation direction. The temperature inhomogeneities 
can be represented by the temperature structure parameter C? 


Neff (1975) and the volume backscatter cross section for 180° 
returns becomes 


2 
oy = 0.0039 Ko (10) 


The acoustic volume scattering cross section is proportional to the 
temperature structure parameter and the acoustic wavenumber. 


This provides the optical structure parameter C® indirectly with high 
spatial and temporal resolution. 
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ll. OBSERVATIONS 


The acoustic sounder was directed vertically into the atmosphere to collect 
atmospheric data from 1 March 2001 until 12 August 2002 at the Starfire Optical 
Range in Kirtland Air Force Base, New Mexico. It was located near the top, but 
on the south side of a 70 m hill. Operating with a 4 KHz, 10-millisecond acoustic 
pulse, the acoustic sounder sampled the temperature structure parameter C? 
between 5 and 150 meters (in height) every second. The optical structure 


parameter C> was later calculated from the measured C* values. The Doppler 


frequencies of the returned signals were also captured. An anemometer located 
15m above the ground measured the wind speeds. 


The raw data from the acoustic sounder and the anemometer were saved 
as data files. Computing the autocorrelation of the mean C,° values between 20- 
50 meters versus time and multiplying by the wind speed gave an autocorrelation 
function versus distance. The acoustic sounder data consisted of C,° profiles 
collected once a second between 20 and 150 meters above the ground with a 
range resolution of 1.7 meters. After computing a single mean C,° value for the 
data between 20-50 meters, the autocorrelation function of these values for 2, 5 
and 10-minute segments of data gave the autocorrelation function versus time. 
The distance between the peak and the first minimum was computed by 
searching for the change in slope, after applying a 5" order, low-pass, 
Butterworth filter to the smooth the irregularities in the autocorrelation function. 
The distance between the peak and the first minimum in the auto correlation 
function, versus time, multiplied by the wind speed gave the thermal 
autocorrelation length. This distance is physically represented as the distance 
between the center of an area of strong thermal activity and the center of the 
area of relatively little thermal activity adjacent to it. Using a Matlab program 
(Appendix A) the data from 1000-1600 (local time) was analyzed. This time 


interval was chosen because this is when convection is intense. The data from 


other times would exhibit different scales and would not represent the strong 


turbulence situations. 


The acoustic sounder data was a continuous, 24-hour stream and had to 
be organized into subsets in order to measure the convective scale lengths. The 
size of the subset chosen affects the measured scale length. For instance, a 
larger time scale would emphasize the distance between individual thermal 
plumes. A time scale of 10 minutes was chosen to measure this length. Smaller 
time scales would reveal features within a thermal plume. Time scales of 2 and 5 
minutes were chosen to look for these features. The Matlab program displays 
this data in a color plot that depicts, in some sense, the plumes themselves. The 
program also uses an autocorrelation function in order to measure the time 
interval between the center of one plume and the center of the quiet region 
between it and the next plume. 


A. RESULTS 


As a representative sample of some of the data, the C> and Doppler 


frequencies measured over a 5-minute or 600-second interval (1 July 2001 from 
1350-1355 hours, MST) are plotted in Figure 3 below. This figure shows a plot of 


the magnitude of C? represented in terms of intensity with the radial Doppler 


frequencies of + 4 m/s as color for a 600 second interval. In general, red regions 
are plumes of air (similar to convective cells) moving upwards, green regions 
have little vertical motion, and blue regions are moving downwards. The 
speckled area mostly above 50 meters in range in Figure 3 represents noise. 
This noise degrades the calculation of the correlation time. Further, the data 
between 0-20 meters was not valid (and is not depicted in the plot because of 
this). For these reasons, the data analyzed were restricted to lie between 20-50 


meters in range. 
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Figure 3. Plot of C> represented as intensity and vertical velocity as color 


between 20 and 150 meters over a 300 second interval (1 July 2001, from 
1350 —1355 MST). 


A technique was devised in the Matlab code to effectively filter out a vast 
majority of the noise in the plots produced. Some noise remained in many 
cases, however, so even with this filtering method being employed, the data was 
still restricted to lie between the limits of 20-50 meters in range. Figure 4 shows 
the same data as Figure 3, but with some of the noise filtered out. Figure 4 also 
shows graphically what the correlation time represents. For this particular case, 
the Matlab program calculated a correlation time of 29 seconds. Only two areas 
are shown in the picture (for clarity) but the autocorrelation function actually 


produced a correlation time based upon the entire 300-second interval. 


Filtered Data 
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Figure 4. Figure 3 with the noise somewhat filtered out. The correlation time 


can be graphically seen as the distance between the center of a plume 
and the center of the quiet, adjacent area. 


Figure 5 shows a plot of the autocorrelation function over the 300-second 
time interval. The Matlab program used this function to determine the correlation 
time. The program found the peak of the curve plotted below, found the first 
minimum to either side of the peak, and found the time between the peak and the 
first minimum. This time is the correlation time and corresponds to the area 


marked off in Figure 4. 


Auto-Correlation of Mean 


0 100 200 300 400 500 600 


Figure 5. Autocorrelation plot of C> mean (1 July 2001, from 1350 -1355 
MST). 


Although it is not so evident in the particular case depicted in Figure 5, 
often the autocorrelation plot would contain too much fine detail and needed to 
be smoothed out in order to find a consistent minimum. A 5" order Butterworth 
filter with a cutoff of W, = 0.2 of the Nyquist frequency was used to accomplish 
this. Figure 6 shows the same plot as Figure 5, but with this Butterworth filter 
applied to it. The plot also depicts the correlation time for this case of 29 


seconds. 


Butterworth Filter on the Correlation of Mean 


Correlation Time: 
29 seconds 





“0 100 200 300 400 500 600 


Figure 6. 5" order Butterworth filter of with W,, = 0.2 applied to the 
autocorrelation plot of C> mean (1 July 2001, from 1350 -1355 MST). 


Depending upon the length of the data sample used to compute a function 
and the strength of the low pass Butterworth filter applied to the autocorrelation 
function, different values for the correlation time were found. Small time samples 
are sensitive to finer features of individual thermal plumes. One can visualize 
this by observing Figure 4 and realizing that as the sample time (the x-axis) 
becomes smaller, the number of distinct plumes in the sample will decrease and 
the individual features of any one plume will become more prominent. 
Conversely, for larger time samples, the correlation time between the plumes, or 
organized clusters, becomes evident. Also, reducing the bandwidth of the low 
pass Butterworth filter, by decreasing W,, smoothes out the autocorrelation 


function thus emphasizing larger scale functions. 


Because the correlation time can vary depending upon the sample time 
and the strength of the Butterworth filter choices, three sample lengths of 2 


minutes, 5 minutes and 10 minutes were analyzed. The Butterworth filter for the 
14 


2 and 5 minute samples was set at W, = 0.2 to maintain full scale structure and 
W, was 0.02 for the 10-minute samples to reveal the major structures. 


Multiplying the correlation time by the wind speed gave the correlation 
distance (the distance between thermal plumes for large time samples, or the 
distance between features of an individual plume for smaller time samples). As 
stated previously, an anemometer provided the wind speed data simultaneously 
with the sounder data. The output from the Matlab program was a list of 
correlation lengths for each day between the hours of 1000 and 1600 MST at 
intervals of 2,5, and 10 minutes. These lists of correlation lengths were then 
statistically analyzed and the results are presented in the following tables and 
figures. 


There are three different standard statistical measures used to analyze a 
collection of data such as these: the mean, the median and the mode. Tables 
1,2 and 3 present the mean, standard deviation, median and mode of the 
correlation distances as found by month for the three different time intervals 


measured. 


Correlation Distance 
2 Minute Summary (W, = 0.2) 


| Mardi] 92 {| 70 | 74 | 50 | 2446 | 
| Aproi] 7s | 63 | 62 | 50 | 2913 | 
| Juloi} 81 | 65 | 60 | 45 | 4522 | 
|Augoi] 73 | 62 | 56 | 35 | 4124 | 
|Sepoi] 7o | 59 | 51 | 35 | 3358 | 
| Octoi] es | 88 | 61 | 38 | 4110 | 
|Novoi] 69 {| 57 | 50 | 35 | 2717 | 
|Decoi] 62 | 67 | 45 | 30 | 3111 | 
| Jano2]| 70 | 60 | 54 | 35 | 3114 | 
| Febo2| 84 | 77 | 62 | 45 | 3320 | 
| Maro2]| 91 | 72 | 69 | 45 | 3750 | 
| Apro2} 99 | 77 | 76 | 68 | 3619 | 
EAE Eee 
Pguto2|—s2—| 6s | ea} a8 | 2640 
|Augo2}| 95 | 97 | 67 | 30 | 779 | 
Bae Ee 


Table 1. Statistical summary of correlation distances for 2-minute interval 
length and W, = 0.2. 
15 





Correlation Distance 
5 Minute Summary (W,, = 0.2) 
| Mar 01) 


| 99 | 76 | 79 | 70 | 985 | 
| Aproi] 84 | 64 | 67 | 50 | 1168 | 
[duno] 72 | 52 | 55 | 45 | 1222 | 
| Juloi] 88 | 59 | 72 | 50 | 1814 | 
|Augoi] 83 | 63 | 67 | 45 | 1650 | 
|Sepoi] 79 | 69 [| 58 | 45 | 1350_| 
| Octoi] 103 | 109 | 69 | 50 | 1649 | 
|Novoi] 80 | 65 | 62 | 45 | 1093 | 
|Decoi] 69 | 65 | 50 | 45 | 1252 | 
| dano2] 80 | 63 | 64 | 45 | 1253 | 
|Febo2] 89 | 71 | 67 | 50 | 1333 | 
|Maro2] 99 | 73 | 76 | 50 | 1486 | 
| Apro2] 101 | 73 | 78 | 53 | 1452 | 
|Mayo2] 101 | 78 | 78 | 45 | 1373 | 
| Juno2z] 92 | 68 [| 71 | 50 | 896 | 
| Julo2] 92 | 67 | 76 | 45 | 1133 | 
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Table 2. Statistical summary of correlation distances for 5-minute interval 
length and W, = 0.2. 





Correlation Distance 
10 Minute Summary (W, = 0.02) 


neo nants] 
| 531 | 


531 498 
p 588 
Jun 0 615 
Jul 0 913 

g 807 
p 679 
Oct 01 796 
553 
| Dec Oi] 374 | 394 | 249 | 188 | 632 
633 
| Febo2| 501 | 406 | 419 | 263 | 671 
| Mar 02] 524 | 369 | 417 | 350 | 756__| 
| Apro2| 592 | 387 | 494 | 375 | 728 | 
|May 02] 552 | 317 | 463 | 450 | 689 | 
| Juno2| 487 | 341 | 391 | 350 | 451 | 
| Aug 02] 473 | 345 | 366 | 303 | 156__| 


a | | eS | eS ee eee 
Entire} io, | 373 301 | 316 | 631 
Period 


Table 3. Statistical summary of correlation distances for 10-minute interval 
length and W, = 0.02. 
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Figure 7 presents histograms of the mean correlation lengths (from Tables 
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Figure 7. Histogram of 2,5, and 10 minute sampled correlation distances for 


March 2001 to August 2002. 


Because of the wide separation in Figure 7 between the peak of the 10- 
minute interval histogram with the other two it is difficult to see where the peaks 
of the two and five minute histograms lie. Figure 8 plots just the two and five 


minute histograms. 
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Figure 8. Histogram of 2-minute and 5-minute sampled correlation distances. 


The statistical mean, median and the mode varied in a consistent manner. 
The shape of the histograms shows why this is so. The histograms all have very 
long tails. This fact skews the median, and especially the mean to the right of the 
most likely value, the mode. The mode is the most likely value to occur and as 
such is deemed the most significant for potential applications of this thesis. 
Figure 9 shows one example of a histogram with the mean, median and mode 
plotted on it. The tail of this histogram is actually much longer than shown - it 
extends to 1140 meters. It has been truncated in this plot in order to better show 


the difference between the mean, median and mode. 
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Figure 9. Histogram of 2-minute sampled correlation distance with mean, 


median and mode plotted. 


Figure 10 displays the mean correlation length, with error bars of plus and 


minus one standard deviation of the data set, for each of the time intervals, by 


month. 
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Mean of the Correlation Distance as a Function of Time of 
Year 
(2-minute intervals) 
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Figure 10. —_ Histogram of 2, 5 and 10-minute mean correlation distances. 
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Figure 11 compares the mean, median and mode for the 2, 5 and 10 
minute intervals. Note that the most likely value, the mode, is also the lowest 


value. 
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Figure 11. | Comparison of the mean, median and mode for the 2, 5 and 10 
minute intervals. 
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Another way to compare the data is examine the means as a group (for 
the three time intervals), the medians as a group and the modes as a group. 
Figure 12 depicts this. 


Means 





2-minutes — — — — 5-minutes <== 10-minutes 














Correlation Distance (m) 
iN 
fo} 
Oo 











s ~ gS SP SH YG YH Y 
w sé we oO X& & oy > SS) 
Date 
Medians 
600 2-minutes ------- 5-minutes 10-minutes 

















Correlation Distance (m) 
wo 
fo) 
pa | 











Modes 











10-minutes 











Correlation Distance (m) 
ny 
oa 
oO 











Figure 12. | Comparison of the means as a group, medians as a group and 
modes as a group for the 2, 5 and 10 minute intervals. 
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It was also of interest to see how the correlation length vaired hourly over 
the course of the day, from 1000 to 1600 MST. The data for the 5-minute interval 
was chosen as a representative sample and the correlation lengths were 
calculated for each hour, by month. Then the months were grouped together into 
seasons: Winter = December, January, February; Spring = March, April, May; 
Summer = June, July, August; and Fall = September, October, November. The 
hourly changes in correlation length were then plotted for each season as well as 


a composite of all the data. This is presented in Figure 13. 
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Figure 13. | Average correlation distance of the seasons, by hour, between the 
hours of 1000 and 1600 MST. 


B. DISCUSSION 
The mean, standard deviation, median and mode were calculated for the 
autocorrelation lengths of the average C,” values for 20-50 meter vertical 


segments of acoustic sounder data grouped in 2, 5 and 10 minute time intervals. 
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These three statistical measures varied from each other because of the long tail 
shape of the histogram illustrated in Figure 9. While each statistical measure 
has value, the mode, which represents the most likely value, is the most 


representative value of the distance between a plume and an adjacent minimum. 


As Figure 12 shows, the mean, median and mode of the 2 and 5 minute 
intervals match up quite well with each other while the 10-minute interval’s 
statistical values are much higher. This is most likely due to two factors. The 
first, and most dominant, reason for the discrepancy is that the bandwidth of the 
low pass Butterworth filter applied to the autocorrelation plot was much narrower 
for the 10-minute case. This was done to suppress the fine scale structure seen 
within a single plume and to reveal the separation between the plumes. This 
longer correlation time, of course, implies a longer correlation distance. The 
second reason that the correlation distance was found to be higher in the 10- 
minute sample was that for longer sample times, larger scale features of the 
thermal plumes were being correlated. A 10-minute segment was long enough 
to include several plumes, while for the 2 and 5-minute intervals, individual 


features within a plume were more prevalent. 


Figure 13, which presents the average correlation lengths of the seasons, 
organized by hour, displays several interesting features. The overall trend of the 
data is that the distance tends to increase from approximately 1000 hours to 
1500 hours and then goes down. In the morning as the sun rises higher in the 
sky, the ground warms up producing thermal plumes. These plumes grow in 
height and lateral extent as they carry more heat up into the boundary layer. The 
length decreases after 1500 hours as the solar radiation abates and the ground 
begins to cool. Spring has the highest correlation distances, which should be 
expected since the average wind speeds in the spring are higher than in the 
other seasons. Fall shows the steepest increase in correlation distance, while 
winter shows a slight decrease. 
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lll. CONCLUSIONS AND RECOMMENDATIONS 


The horizontal separation between convective thermal plumes, and 
features within a thermal plume, were determined from acoustic sounder profiles 
with the corresponding anemometer wind speeds and extensive data analysis. 
The mean, standard deviation, median and mode were calculated for the 


computed correlation lengths. The acoustic sounder produced a vertical C* 
profile once a second that was represented as a single C* average for the 20- 


50m altitude range. These averages were grouped into 2, 5 and 10 minute long 
sequences to compute an autocorrelation function. The times between the 
autocorrelation maxima and the first minima were found and multiplied by the 
average wind speed over that interval. The data sampled at 2 and 5-minute 
intervals most likely revealed features within an individual thermal plume. The 2 
and 5-minute time interval statistics correlated very strongly with each other, and 


the mean correlation distance was 81 meters + 70 meters and 89 meters + 72 


meters, respectively. Their medians were 61 meters and 69 meters; and their 
modes were 41 meters and 50 meters, respectively. These two time intervals 
also used the same strength Butterworth low pass filter in the Matlab program to 
determine the correlation time, which provided just enough smoothing to find a 


robust minimum in the autocorrelation function. 


The 10-minute time interval statistics were significantly larger than the 2 


and 5-minute cases. The mean correlation length was 494 meters + 373 meters, 


the median was 391 meters and the mode was 316 meters. This was the result 
of using a narrower bandwidth Butterworth filter to smooth the autocorrelation 
function. This stronger filter was designed to filter out the smaller features within 
an individual plume and, therefore, find the distance between the center of a 
plume and the center of the quiet region adjacent to that plume (effectively 
finding the distance between thermal plumes). Another factor involved in the 


higher values for the 10-minute interval was the fact that the time interval was 
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longer: 10 minutes versus 2 or 5 minutes. This longer time interval allowed 


several thermal plumes to be correlated in the same data set. 


The validity of this study should be verified using optical methods to 
determine the correlation lengths between and within thermal plumes. This 
would be important to do before any major applications utilizing the results 
obtained here were begun. If the results agree, then this acoustic sounder 
method could be a cost effective way to study thermal characteristics, and the 


resulting degradation in the optical environment, at other locations. 
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APPENDIX MATLAB CODE 


The following program is the code used to generate the data analyzed in 


this thesis. 


% This program reads in acoustic sounder data and 
% meteorological data collected at the Starfire Optical 
% Range (SOR) and produces an output listing in successive columns: 


% year of sample 

% month of sample 

% day of sample 

% time of sample 

% correlation time (unfiltered), 

% correlation time (filtered) 

% average windspeed (for time of sample) 
% correlation distance (max to first min) 
% correlation distance (max to global min) 


% The video display mode must be set for 24 bits or more 
% 25 Feb 98; 24 Oct 97, DLW 


OEE OES E SERIE EIEIO IO IO AO EO TAIT IIR III IIIT I 
0 


% This file modified Aug/Sep 2004 by Jeff Douds * 


C aalaalaalaaalaelaailabebaaeebaebaeiaeieliiaelibaelebabelabaebaiae 
) 


OEE HERESIES IE IEEE AO IO AO EET TI TID TIO TIO III III IIIT TIT TIT TIT TID AI AIA 
0 





% SOR VARIABLES ON DATA FILES (included for reference) ‘ 
% fdat(1)=floor(now); % days is 
% fdat(2)=rem(now,1); % fraction of day i 
% fdat(3)=Pres; a 
% fdat(4)=Temp; s 
% fdat(5)=Rel; ; 
% fdat(6)=atten2; - 
% fdat(7)=Ct_cal; . 
% fdat(8)=Pltrng; * 
% fdat(9)=dt; ‘i 
% fdat(10)=Cn_sensor; e 
% fdat(11:Pltrng+10)= CT%2 array in K*2 m’-2/3 ™ 
% fdat(Pltrng+11:2*Plitrng+10)= Radial Doppler wind speed in m/s . 
% 
% METEOROLOGICAL DATA VARIABLES (in txt files) = 
% col 1=ABD temp (deg C) , 
% col 2 = ABD dew point temp (deg C) * 
% col 3 =ABD wind direction (deg) ig 
% col 4 =ABD wind speed (mph) « 
% col 5 =ABD pressure (mb) . 
% col 6 =ABD pyranmometer Z 
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% 
% 
% 
% 
% 


col 7 = 3.5m wind speed (mph) 


col 8 = 3.5m temp (deg C) 


col 9 = 3.5m dew point temp (deg C) 


col 10 = 3.5m wind direction (deg) 


col 11 = time (hhmmss.00) (increases in 1 minute increments) 


CA atalatalaelaebalaeiaebiaelaeebaeebaebaeiebiaebiaebiaeebeebebelielebaldelabaebehaebebaetehaetebaitebel 
0 


clear all; 


close all; 


cle; 


% ‘0 


Set Constants 


color_scale=1/3/4; 
Ct_cor=1; 
lambda_o=0.5e-6; 
romin=15; 
ro_Pltrng=200; 
scale=2000; 


% ‘0 


User Inputs 


% (1/3)/5 For full scale of 5 m/ 

% Empirical constant (use 2.3 before Feb 10 98) 
% Optical wavelength 

% ro minimum array index altitude (ie 5 = 4m) 

% ro plot maximum range 


% colormap CT“2 scale factor 


time_samp=input(‘\nEnter sampling time in min = '); 


nscans= time_samp*60;% max # of pulses displayed on plot 


tstart=input('Enter time to start (use UTC, in whole hours, e.g. 17)'); 


tend=input('Enter time to end (use UTC, in whole hours, e.g. 23) '); 


% Derived Constants 


nscan1=nscans-1; 
ro_const=1000*2.1*(1.46*(2*pi/lambda_o)*2)4-0.6; 


for j=1:12 %select the month (run through them all sequentially) 


switch j 
case 1 
month='Jan'; 
case 2 
month='Feb'; 
case 3 
month='Mar'; 
case 4 
month='Apr'; 
case 5 
month='May'; 
case 6 
month='Jun'; 
case 7 
month='Jul'; 
case 8 
month='Aug'; 
case 9 
month='Sep'; 
case 10 
month='Oct'; 
case 11 
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* 


month='Nov'; 
case 12 
month='Dec'; 
end 


% loop to read in a month at a time day by day 
for count=1:31 
if count<10 
count=(['0',int2str(count)]); 
else 
count=([int2str(count)]); 
end 


%assign names to the input files 
sounder_file=(['c:\SOR FOLDER\SOR DATA FILES\',month,' 2002\',count,'-', month, '-2002.bin']); 
met_file1=(['c\SOR FOLDER\S2002\',month,'02\'",count,month,'02data.txt']); 


%open the input files 
fid=fopen(sounder_file,'r'); 
met_file=load(met_file1); 


%open the output file 
outputfile=(['c:\SOR FOLDER\Output\'",month,'_2002_10_min_intervals_no_douds_filter.txt']); 
fid2 = fopen(outputfile,'a'); 





%check for valid data files 

ValidData=1; % 1 = yes, valid data; 0 = no, not valid data 

if size(met_file)==[0,0] | fid==-1 
fprintf(‘The file %s does not exist.\n',sounder_file); %output to data file 
ValidData=0; %flag for invalid data 

end 


if ValidData== %only enter this loop if two valid data files have been read 


% read input parameters 
fdat=zeros([10 1]); 
fdat=fread(fid, 10,'float32'); 
Pltrng=fdat(8);  % get vector length 
fsize=2*Pltrng+10; % actual file size 
frewind(fid); 
fdat=zeros([fsize 1]); % Output file data buffer 
fdat=fread(fid,fsize, ‘float32'); 


% Set Arrays 
i0=0:Pltrng-1; 
i1=1:Pltrng; 
n2=1:nscans; 
avg=zeros(size(i1)); 
Cn2=avg'; 
Cn2_row=zeros(Pltrng,nscans, 1); 
Cn2_norm=zeros(Pltrng,nscans, 1); 
color=avg; 
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img=uint8(zeros(Pltrng,nscans,3)); % plot array with 1m range bins 
img2=uint8(zeros(Pltrng,nscans,3)); % plot array with 1m range bins (for 2nd color plot) 

inten=avg'; 

ravg=avg'; % Cn*2 running average 

rgb=zeros(Pltrng,3); % Current data array 
rgb2=zeros(Pltrng,3); % for second color plot 

vavg=avg; 

vdop=avg; 

roi=zeros(1,nscans); 


% Initialize first plot 
close(‘all','hidden’) 
figure('Position’,[430 50 400 300)); 
colormap(hsv); 
m=colormap; 
h1=image(img,'Erasemode','none’); 
axis xy; % set plot for Cartesian coordinates 
xlabel('Time (sec)'); 
ylabel('Range (m)'); 
title(‘Filtered Data’); 
YLim([20 150]); %don't display the bottom 20 meters since this data is not being used anyway 


%lnitialize second plot 
figure(2) 
set(2,'Position’,[20 50 400 300)); 
colormap(hsv); 
m2=colormap; 
h2=image(img,'Erasemode','none’); 
axis xy; % set plot for Cartesian coordinates 
xlabel('Time (sec)'); 
ylabel('Range (m)'); 
title(‘Unfiltered Data’) 
YLim([20 150]); %don't display the bottom 20 meters since this data is not being used anyway 


% Queing loop: this loop queues the data to begin reading at the tstart input 
fdat(2)=0; 
if size(fdat)==[fsize 1] %only enter the following loop if data has been read each time 
while ValidData==1 & (tstart>24*fdat(2)), Y%flag is to check that data has actually been read each time 
fdat=fread(fid,fsize,'float32'); 
if size(fdat)==[fsize 1] 


else 
ValidData=0; 
end 
end 
end 


% Main Data Acquisition Loop (goes until tend is reached) 
for time=tstart:time_samp/60:tend, 


% Loop for each frame that is printed out 
for k=1:nscans 
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fdat=fread(fid,fsize, ‘float32'); 


ValidData==1 & size(fdat)==[fsize 1] %only enter the following loop if data has been read each time 
Pres=fdat(3); 
Temp=fdat(4)+273.15; 
Cn_cal=Ct_cor*(79e-6*Pres/Temp’2)*2; 

avg=(fdat(11:Pltrng+10)); 

inten=min(scale*sqrt(max(avg,0)),255); 

inten_unfiltered=inten; 


%the following for loop will black out the areas above where the avg 
%Vvalue was first <=0 and set the avg values above there to zero 
%the avg values at these points are already very close to 
%or less than zero 
flag=0; 
for counter=1:Pltrng 
if or((avg(counter)<0), (flag==1)) 
inten(counter)=0; 
avg(counter)=0; 
flag=1; 
end 
end 


%the following algorithm will look for noise beginning in the upper 
%region and working its way down. As it finds noise, it will set avg 
Y%eand inten to zero 
noise=1; 
for j=150:-1:2 %look from the to of the column down to the bottom 
if noise==1; 
if avg(j)>0.002 %avg value too high ==> noise; set to zero 
avg(j)=0; 
inten(j)=0; 
elseif avg(j)==0 %%avg value=0 is inconclusive, so check the next pt down the column 


else 
noise=0; %once avg value drops below threshold level, back out of the loop 
%and stop changing avg and inten values. This will preserve the valid 
%high avg numbers at the bottom of the columns. 
end 
end 
end 


% Doppler processing 
vdop=transpose(fdat(Pltrng+1 1:2*Pltrng+10)); 
color=max(min(color_scale*vdop+0.333,0.666),0); % Clip red & blue 
color=round(color*63+1); % Color table index 
rgb(i1,:)=m(color(i1),:); 
for l=1:3,rgb(:,l)=inten.*rgb(:,l);end 


% Shift plot and display 
img(:,1:nscans-1,:)=img(:,2:nscans,:); 
img(:,nscans, :)=uint8(round(rgb)); 
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ravg=0.6*ravg+0.4*avg; 

Cn2=Cn_cal*avg; 
Cn2_row(:,1:nscans-1)=Cn2_row(:,2:nscans); 
Cn2_row(:,nscans)=Cn2(:); 
row_avg=mean(Cn2_row,2); 
for i_col=1:nscans, 

Cn2_norm(:,i_col)=Cn2_row(:,i_col)-row_avg(:); 





end 


%Produce data for 2nd color plot 
rgb2(it ,:)=m2(color(i1),:); 
for l=1:3,rgb2(:,l)=inten_unfiltered.*rgb2(:,l);end 


% Shift plot and display 
img2(:,1:nscans-1,:)=img2(:,2:nscans,:); 
img2(:,nscans, :)=uint8(round(rgb2)); 


end %if loop to check if fdat had data in it 
end % loop for each frame that is printed out 
if size(fdat)==[fsize 1] %only enter the following inner loop if data has been read each time 


% prepare the output for the color graph and draw it 
set(h1,'CData',img) 
drawnow 


% prepare the output for the color graph and draw it 
set(h2,'CData',img2) 
drawnow 


% find the correlation time between limits in array below 
Cn2_ab = mean(Cn2_norm(20:50,:)); %Find mean of Cn2 (finds the mean value of each column) 
cor_ht1 = xcorr(Cn2_ab); %Find correlation time 


% Filter the data 
[b,a] = butter(5,0.02); %5th order, Wn=0.2 
corfilt = filtfilt(b,a,cor_ht1); 


% Check for first min to max 
Ne=nscans; 
last=cor_ht1(nscans); 
while ((corfilt(N-1)<last) & (N>2)), 
last=corfilt(N-1); 
N=N-1; 
end 


time_O=nscans-N; 
% Check for global min 


[¥1,11] = min(corfilt(1 :nscans)); 


time_1 = nscans-l1; 
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% Find the corresponding time in the met data file 
met_tstart_1=time*10000; %decimal format (e.g. 17.8333) 
hours=double(int8(met_tstart_1/10000)); 
met_tstart_2=met_tstart_1 - hours*10000; %leaves only the minutes and seconds, still in base 10 format 
minutes=double(int8(met_tstart_2/100)); %extracts the minutes 
minutes=double(int8(minutes*0.6)); convert minutes to time format 
seconds=double(met_tstart_2 - (double(int8(met_tstart_2/100))*100)); %extracts the seconds and converts to time format 
seconds=double(int8(seconds*0.6)); %Yconvert seconds to time format 
met_tstart=(hours*10000) + (minutes*100) + seconds; %time in the correct format for the .txt met file 
J=min(find(met_file(:,11)>=met_tstart)); %finds the row in the met file where the time starts 
sum_wndspd=0; 
for i=J:J+(time_samp — 1) 
sum_wndspd=sum_wndspd + met_file(i,7); Yadd the windspeeds over the sample to then calculate an average 
end 
avg_wndspd=sum_wndspd/time_samp; “using the 3.5m windspeed, not the ABD windspeed 


% calculate the correlation distance 
cor_dist_first=0.44704*avg_wndspd*time_0; %0.44704 = conversion factor for mph to m/s; distance from max to first min 
cor_dist_global=0.44704*avg_wndspd*time_1; %distance from max to global min 


figure(3); 

set(3,'Position’,[20 430 400 320)); 
plot(cor_ht1); 
title(‘Auto-Correlation of Mean’); 


figure(4); 

set(4,'Position’,[430 430 400 320)); 

plot(corfilt); 

title(‘Butterworth Filter on the Correlation of Mean’); 


Tm=datevec(fdat(1)+fdat(2)); % pulls out the year, month, day from the SOR data file 
% output to screen 
fprintf(‘\nYear\tMonth\tDay\t\tUTC Time\tCorr Time 1\t\tCorr Time 2\t\tAvg Wndspd\t\tCorr Dist 1\t\tCorr Dist 2'); 
fprintf(‘\n%g\t Yg\t\t Yg\t\t %02.0F:%02.0F:%02.0f\t %4.Of\t\t\t %4.Of\t\t\t %8.2f\t\t %8.2f\t\t %8.2f\t\r\n’.... 
Tm(1),Tm(2),Tm(3),hours,minutes,seconds,time_0,time_1,avg_wndspd,cor_dist_first,cor_dist_global); 

if Tm(1)==2001 | Tm(1)==2002 %some bad data has the incorrect date, this will filter that data out by not printing it 

% output to data file 

fprintf(fid2,'%g\t Yg\t Yg\t %8.4f\t %02.0F:%02.0f:%02.0f\t %4.0f\t %4.0f\t %8.2f\t %8.2f\t %8.2fi\t\r\n’,... 

Tm(1), Tm(2),Tm(3),time,hours,minutes,seconds,time_0,time_1,avg_wndspd,cor_dist_first,cor_dist_global); 
end 
end %inner if loop to check if fdat had data in it 
pause; 
end % Main Data Acquisition Loop 


fclose(fid); % close current SOR file 


end %if statement for valid data obtained 
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end % for loop to read in a day at time for the whole month 


end %for loop to select the month 


fclose(fid2); % Close output file 
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